Nonlinear damping in a micromechanical oscillator 
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Nonlinear elastic effects play an important role in the dynamics of microelectromechanical systems 
(MEMS). A Duffing oscillator is widely used as an archetypical model of mechanical resonators with 
nonlinear elastic behavior. In contrast, nonlinear dissipation effects in micromechanical oscillators 
are often overlooked. In this work, we consider a doubly clamped micromechanical beam oscillator, 
which exhibits nonlinearity in both elastic and dissipative properties. The dynamics of the oscillator 
is measured in both frequency and time domains and compared to theoretical predictions based on 
a Duffing-like model with nonlinear dissipation. We especially focus on the behavior of the system 
near bifurcation points. The results show that nonlinear dissipation can have a significant impact 
on the dynamics of micromechanical systems. To account for the results, we have developed a 
continuous model of a geometrically nonlinear beam-string with a linear Voigt-Kelvin viscoelastic 
constitutive law, which shows a relation between linear and nonlinear damping. However, the exper- 
imental results suggest that this model alone cannot fully account for all the experimentally observed 
nonlinear dissipation, and that additional nonlinear dissipative processes exist in our devices. 



I. INTRODUCTION 

The field of micro-machining is forcing a profound re- 
definition of the nature and attributes of electronic de- 
vices. This technology allows fabrication of a variety of 
on-chip fully integrated micromechanical sensors and ac- 
tuators with a rapidly growing range of applications. In 
many cases, it is highly desirable to shrink the size of 
mechanical elements down to the nano-scale jl|, [2, H, 
This allows enhancing the speed of operation by increas- 
ing the frequencies of mechanical resonances and improv- 
ing their sensitivity as sensors. Furthermore, as devices 
become smaller, their power consumption decreases and 
the cost of mass fabrication can be significantly lowered. 
Some key applications of microelectromechanical systems 
(MEMS) technology include magnetic resonance force 
microscopy (MRFM) [5, 6] and mass-sensing S, & • 
Further miniaturization is also motivated by the quest 
for mesoscopic quantum effects in mechanical systems 

[m, E H m E H E [3. 

Nonlinear effects are of great importance for microme- 
chanical devices. The relatively small applied forces 
needed for driving a micromechanical oscillator into the 
nonlinear regime are usually easily accessible Thus, 
a variety of useful applications such as frequency synchro- 
nization [20l|, frequency mixing and conversion (2ll . [22| . 
parametric and intermodulation amplification [23|, me- 
chanical noise squeezing [24], stochastic resonance [25*], 
and enhanced sensitivity mass detection (26| can be im- 
plemented by applying modest driving forces. Further- 
more, monitoring the displacement of a micromechanical 
resonator oscillating in the linear regime may be difficult 
when a displacement detector with high sensitivity is not 



available. Thus, in many cases the nonlinear regime is 
the only useful regime of operation. 

Another key property of systems based on mechani- 
cal oscillators is the rate of damping. For example, in 
many cases the sensitivit y o f NEMS sensors is limited 
by thermal fluctuation |l,l23|, which is related to damp- 
ing via the fluctuation dissipation theorem. In general, 
micromechanical systems suffer from low quality factors 
Q relative to their macroscopic counterparts [1, [H, [2§|. 
However, very little is currently known about the un- 
derlying physical mechanisms contributing to damping 
in these devices. A variety of different physical mech- 
anisms can contribute to damping, including bulk and 
surface defects [sO*, "sT], thermoelastic damping (32I. [33|. 
nonlinear coupling to other modes, phonon-electron cou- 
pling, clamping loss [H, [35|, interaction with two level 
systems [36], etc. Identifying experimentally the con- 
tributing mechanisms in a given system can be highly 
challenging, as the dependence on a variet y o f parame- 
ters has to be examined systematically (stI. Issl. Isol. [4Q|. 

The archetypical model used to describe nonlinear 
micro- and nanomechanical oscillators is the Duffing os- 



cillator 41 1 . This model has been studied in great depth 
lijiill, and special emphasis has been given to 
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the dynamics of the system near the bifurcation points 

0, El, 53, S, SlHli. 

In order to describe dissipation processes, a linear 
damping model is usually employed, either as a phe- 
nomenological ansatz, or in the form of linear coupling to 
thermal bath, which represents the environment. How- 
ever, nonlinear damping is known to be significant at 
least in some cases. For example, the effect of nonlinear 
damping for the case of strictly dissipative force, being 
proportional to the velocity to the n'th power, on the re- 
sponse and bifurcations of driven Duffing fHll. [52I. IHsl . [53 | 
and other types of nonlinear oscillators fSl, [53, [55| has 
been studied extensively. Also, nonlinear damping plays 
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an important role in parametrically excited mechanical 
where without it, solutions will grow with- 
57 

In spite of the fact that a massive body of litera- 
ture exists which discusses the nonlinear elastic effects 
in micro- and nanomechanical oscillators as well as the 
consequences of nonlinear damping, the quantitative ex- 
perimental data on systems with nonlinear damping, es- 
pecially those nearing bifurcation points, remains scarce. 
Furthermore, such systems impose special requirements 
on the experiment parameters and procedures, mainly 
due to the very slow response times near the bifurcation 
points. Straightforward evaluation of these requirements 
by simple measurements can facilitate accurate data ac- 
quisition and interpretation. 

In the present paper we study damping in a microme- 
chanical oscillator operating in the nonlinear regime ex- 
cited by an external periodic force at frequencies close 
to the mechanical fundamental mode. We consider a 
Duffing oscillator nonlinearly coupled to a thermal bath. 
This coupling results in a nonlinear damping force pro- 
portional to the velocity multiplied by the displacement 
squared. As will be shown below, this approach is equiv- 
alent to the case where the damping nonlinearity is pro- 
portional to the velocity cubed [58]. In conjunction with 
a linear dissipation term, it has also been shown to de- 
scribe an effective quadratic drag term [59]. 

We find that nonlinear damping in our micromechani- 
cal oscillators is non-negligible, and has a significant im- 
pact on the oscillators' response. Furthermore, we de- 
velop a theoretical one-dimensional model of the oscilla- 
tor's behavior near the bifurcation point [H, [45|. Most 
of the parameters that govern this behavior can be es- 
timated straightforwardly from frequency response mea- 
surements alone, not requiring exact measurement of os- 
cillation amplitudes. Measuring these parameters under 
varying conditions provides important insights into the 
underlying physical mechanisms [60, £1] . 

We use our results to estimate different dynamic pa- 
rameters of an experimentally measured micromechanical 
beam response, and show how these estimations can be 
used to increase the accuracy of experimental measure- 
ments and to estimate measurement errors. The main 
source of error is found to be the slowing down behav- 
ior near the bifurcation point, also known as the saddle 
node "ghost" [46.]- We also investigate the possibility of 
thermal escape of the system from a stable node close to 
the bifurcation point [ll, [13, S, ^ and find that the 
probability of this event in our experiments is negligible. 

Finally, we propose and analyze a continuum mechan- 
ics model of our micromechanical oscillator as a pla- 
nar, weakly nonlinear strongly pretensioned, viscoelas- 
tic beam-string. The analysis of this model illustrates 
a possible cause for non negligible nonlinear damping as 
observed in the experiment. 



II. EXPERIMENTAL SETUP 

For the experiments we employ micromechanical os- 
cillators in the form of doubly clamped beams made of 
PdAu (see Fig. [T]). The device is fabricated on a rectan- 
gular silicon-nitride membrane (side length 100-200 pm) 
by the means of electron beam lithography followed by 
thermal metal evaporation. The membrane is then re- 
moved by electron cyclotron resonance (ECR) plasma 
etching, leaving the doubly clamped beam freely sus- 
pended. The bulk micro-machining process used for sam- 
ple fabrication is similar to the one described in The 
dimensions of the beams are: length 100-200 pm, width 
0.25-1 pm and thickness 0.2 pm, and the gap separating 
the beam and the electrode is 5-8 pm. 

Measurements of all mechanical properties are done 
in-situ by a scanning electron microscope (SEM) (work- 
ing pressure 10~^ Torr), where the imaging system of the 
microscope is employed for displacement detection |ll6| . 
Some of the samples were also measured using an op- 
tical displacement detection system described elsewhere 
(23 |. Driving force is applied to the beam by applying a 
voltage to the nearby electrode. With a relatively mod- 
est driving force, the system is driven into the region of 
nonlinear oscillations [16, 63]. 




FIG. 1. A typical device consists of a suspended doubly 
clamped narrow beam (length 200 pm, width 1-0.25 pm, and 
thickness 0.2 pm) and a wide electrode. The excitation force is 
applied as voltage between the beam and the electrode, (a) 
Experimental setup and typical sample's dimensions. The 
direction of the vibration of the micromechanical beam is de- 
noted by dotted arrow, (b) SEM micrograph of a device with 
one wide electrode and two narrow doubly clamped beams. 

We use a network analyzer for frequency domain mea- 
surements, as shown in Fig.[2l For time domain measure- 
ments of the slow varying envelope we employ a lock-in 
amplifier, connected as show in Fig. [3l The mechani- 
cal oscillator is excited by a monochromatic wave, whose 
amplitude is modulated by a square wave with low fre- 
quency (20-50 Hz). This results in bursting excitation, 
which allows measurement of ring-down behavior in time 
domain. The lock-in amplifier is locked to the excitation 
frequency, and measures the amplitude of the slow en- 
velope of the oscillator's response. The lock-in amplifier 
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time constant should be much smaller than the ring down 
time, which is governed by dissipation in the microme- 
chanical system. Typically, in our experiments, the time 
constant is 100 ps and the characteristic ring down time 
is 10 ms. 
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FIG. 2. Network analyzer is used for frequency domain mea- 
surements. If the system is excited into a bistable regime, 
special care should be taken to ensure accurate measurement 
near bifurcation points, as discussed in Sec. IIIIEI 



Electron Detector 
beam 



Lock In Amplifier 



Signal 
In 




DC + 



FIG. 3. Lock-in amplifier is employed for time domain mea- 
surements. The oscillator is excited at a single frequency. The 
amplitude of the excitation is modulated by a square wave, 
effectively turning the excitation on and off 20-50 times per 
second. Such bursting excitation is used to measure the ring- 
ing down of the slow envelope in the time domain. 

The displacement detection scheme described above is 
not exactly linear, because the amount of the detected 
secondary electrons or reflected light is not strictly pro- 
portional to the mechanical oscillator amplitude, but 
merely a monotonic function of the latter. Nonuniform 
distribution of primary electrons or light power in the 
spot increases this nonlinearity even further. Thus, some 
distortion in the measured response amplitude is intro- 
duced. 



III. THEORY 
A. Equation of motion 

We excite the system close to its fundamental mode. 
Ignoring all higher modes allows us to describe the dy- 
namics using a single degree of freedom x. 

In the main part of this study, no assumptions are 
made about the source of linear and nonlinear dissipa- 
tion. The energy dissipation is modeled phenomenolog- 
ically by coupling the micromechanical oscillator to a 
thermal bath consisting of harmonic oscillators [HI, [H, 
IggI . [67| . Physically, several processes may be responsi- 
ble for mechanical dampi ng l27.. .281. Is^ . IgsI . [6§| , includ- 
ing thermoelastic effects [32|, [SJ, [70|, friction at g rain 
boundaries [71], bulk and surface impurities (29l. It^TItsI. 
electrical losses, clamping loss [35, tI, [zll, etc. We also 
regard the linear and nonlinear damping constants as in- 
dependent of one another, although they probably result 
from same physical processes. In Sec. IV Bl we consider 
one possible model connecting the linear and nonlinear 
dissipation coefficients, and compare its predictions to 
experimental data. 

The Hamiltonian of the system, which includes the me- 
chanical beam and thermal bath modes coupled to it, is 



(1) 



where 



2m 



V 2m5 



Hi = y^^T{x,Ub)qb, 



describe the micromechanical beam, the thermal bath, 
and the interaction between them, respectively. Here, m 
is the effective mass of the fundamental mode of the mi- 
cromechanical beam, and p and x are the effective mo- 
mentum and displacement of the beam. Also, U{x) is 
the elastic potential, and Scap{x,t) = C{x)V{t)'^ /2 is the 
capacitive energy, where C{x) = Co / {1 — x / d) is the dis- 
placement dependent capacitance, d is the gap between 
the electrode and the beam, and V{t) is the time de- 
pendent voltage applied between the electrode and the 
micromechanical beam. The sum denotes summing 
over all relevant thermal bath modes, while cOb is the fre- 
quency of one of the modes in the thermal bath with 
effective momentum pij and displacement ^5, and 7715 is 
the effective mass of the same mode. Finally, r(x,cj5) 
is a function describing the interaction strength of each 
thermal bath mode with the fundamental mode of the 
micromechanical beam. 
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The equations of motion resulting from ([T]) are 



Qb- 



dx 



The formal solution of ([2bl) can be written as 



(2a) 
(2b) 



Qb\t) = Qbo cosiJbt H smcjbt 

UJb 



^ T{x,UOb]T) 

rribUJb 



sina;5(r — t)(ir, 



or, integrating by parts, 



^6(^) = ^60 COS cobt + — sin cobt 

UJb 



r(x,cj5;0) ^ r{x,ujb;t) 

2 COS Ubt 



rribUJ^ 

' x{r) dT{x,u;b]r) 
mbool dx 



rubujl 



coscobir — t)dr^ (3) 



where Qbo and qbo are the initial conditions of the 
thermal mode displacement and velocity, respectively; 
and r{x^ujb;s) denotes the coupling strength function 
r{x^ujb) evaluated at time s. 

Substituting (|3]) into ([2a|) . one gets 



mx -\- / /C(x, r)x(r)(i7 
Jo 



where n{t) is the noise. 



dU{x) 
dx 

d£caip{x^ i) 

dx 



mn(t), (4) 



U{x) = U{x) - ^ 



2mbu;? 



(5) 



is the renormalized potential, and 

9r(x, UJb] t) 9r(x, UJb; t) cosuJb{r — t) 



JC{x,t,r) = 



dx 



dx 



rribujl 



is the memory kernel [13, Also, the initial slip term 
given by 

Q 

^ r(x, UJb] 0) cosuJbt—T{x, UJb; t)/ (rribUjl), 



has been dropped [67|. Finally, the noise autocorrelation 
for an initial thermal ensemble is 

{n{t)n{s)) = ■^^/C(x,t,5), 

where T is the effective temperature of the bath, and 
/cb is the Boltzmann's constant. The last result is a 



particular form of the fluctuation-dissipation theorem 

H [ri, [7i,[83|. 

We employ a nonlinear, quartic potential U{x) = 
\kix'^ + jksx^ in order to describe the elastic proper- 
ties of the micromechanical beam oscillator. Assuming 
r(x, UJb) to be polynomial in x, it can be deduced from ([5|) 
that only linear and quadratic terms in T{x^uJb) should 
be taken into account [66, 8l|, i.e., 

T{x,ujb) = giiujb)x ^ ^g2{ujb)x'^ . (6) 
The memory kernel in this case is 

/C(x, t,r) (^gf + gig2{x{t) + x{t)) 

b 

^glx{t)x{r) 



rribUjl 



Making the usual Markovian (short-time noise autocor- 
tion m, m 

S{t — s), one obtains 



relation) approximatior 



I, i.e., /C(x,t,5) (X 



JC{x, t, r) = (26ii + b2X + bsix^) S{t - r), 
and the equation of motion dU becomes 



mx + (26ii + 631^^ + bs2x'^)x + kix + ksx^ 



dx 



mn(t), (7) 



where 611 is the linear damping constant, 631 and 632 are 
the nonlinear damping constants, ki is the linear spring 
constant and k^, is the nonlinear spring constant. 

Some clariflcations regarding ([7j) are in order. The 
quadratic dissipation term h2xx has been dropped from 
the equation because it has no impact on the flrst or- 
der multiple scales analysis, which will be applied below. 
An additional dissipation term proportional to the cubed 
velocity, 632^^, has been added artiflcially. Such term, 
although not easily derived using the analysis sketched 
above, may be required to describe some macroscopic 
friction mechanisms [4l|, [5l|, [55|, such as losses associ- 
ated with nonlinear electrical circuits. It will be shown 
below that the impact of this term on the behavior of the 
system is very similar to the impact of bsix'^x. 

The applied voltage is composed of large constant (DC) 
and small monochromatic components, namely, V{t) = 
^DC -\- V cos ujt. The one dimensional equation of motion 
(0 can be rewritten as 



X + (2711 + 731X + 732X )x + ujqX + asx 

^DC ' 2 



Co {V§c + h^'^ ^ 2Vbcv COS ujt^ cos2ujt) 



2md{l-^J 



n{t), (8) 



where cjg = ki/m, 711 bu/m, 731 631/m, 732 
^32/^, and as = ks/m. 
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B. Slow envelope approximation 

In order to investigate the dynamics described by the 
equation of motion (|8]) analyticahy, we use the fact that 
nonhnearities of the micromechanical oscihator and the 
general energy dissipation rate are usuahy smah (as 
shown in Sec. HVl the hnear quahty factor in our sys- 
tems has a typical value of several thousands). In the 
spirit of the standard multiple scales method [4l|, (H], 
we introduce a dimensionless small parameter e in (|8]), 
and regard the linear damping coefficient 711 = 6711, the 
nonlinear damping coefficients 731 = 6731 and 732 = 6732, 
the nonlinear spring constant as = eas, and the excita- 
tion amplitude v = ev diS small. It is also assumed that 
the maximal amplitude of mechanical vibrations is small 
compared to the gap between the electrode and the me- 
chanical beam d, i.e., x/d = ex/d. Also, the frequency of 
excitation co is tuned close to the fundamental mode of 
mechanical vibrations, namely, = cc^o + cr, where a = ea 
is a small detuning parameter. 

Retaining terms up to first order in e in (|8]) gives 



where a is a complex amplitude and c. c. denotes complex 
conjugate. The "slow varying" amplitude a varies on a 
time scale of order Ti or slower. 

The secular equation [4l|, (H] , which follows from sub- 
stitution of (dH) into (jlObp . is 

2uo {jd + (j7i - Au) a) + (3^3 + jjs^o) a^a"" = /o, 



where 



71 = 711 + 731 



and 



where 



73 = 731 + 3cJo732, 



F f 3F 1 



(12) 

(13a) 
(13b) 



(14) 



X + UOqX 

+ e (2711 +73ix^ +732i:^)x + a3X^ - 

= F + 2e/ocoscjt, (9) 



— x{x + ^qX) 

d 



where F = C^V^fj / {2md), and e/o = /o = 
C'oVDC^/(2m(i). We have dropped the noise from the 
equation of motion, and will reintroduce its averaged 
counterpart later in the evolution equation (p!5]) . 

Following [44], we introduce two time scales Tq = t and 
Ti = et, and assume the following form for the solution: 

=xo(To,Ti)+exi(To,Ti). 

It follows to the first order in e that 

d 



d 
di 



d 



and di]) can be separated according to different orders of 



e, giving 



and 



d'^xo 



- CJqXo 



(10a) 



d'^xi 



■ lJqXi = 2/0 cos(cjo^o + ^Ti) 



o~ , ~ 2 , ~ fdxo\ \ dxo .3 

2711 + 731^0 + 732 ^ ] ^^^0 



dToJ J dTo 



represents a constant shift in linear resonance frequency 
due to the constant electrostatic force F. Equation ([T2|) 
is also known as evolution equation. Note that we have 
returned to the full physical quantities, i.e, dropped the 
tildes, for convenience. Also, one must always bear in 
mind that the accuracy of the evolution equation is lim- 
ited to the assumptions considered at the beginning of 
this Section. 

As was mentioned earlier, both nonlinear dissipation 
terms give rise to identical terms in the evolution equa- 
tion (p!2]) . Therefore, the behavior of these two dissi- 
pation cases is similar near the fundamental resonance 
frequency ujq- Also, note that linear dissipation coeffi- 
cient 71 (|13ap is not constant, but is rather quadratically 
dependent on the constant electrostatic force F due to 
the nonlinear dissipation term 731. 

The secular equation (p!2|) can be written as 

jd + (j7i - Au)a + g(l + jp)a^a* = (/o + nsiow(^)) , 

(15) 

where dot denotes differentiation with respect to (slow) 
time, 



Q - 
p- 



3^3 

2cjo' 
73^0 

3^3 ' 



(16) 
(17) 



and nsiow(^) is the av erag ed noise process with the fol- 
lowing characteristics [42|, [43| : 



(^slow(^)) 0, 

(^slow(^)^slow(5)) = N5{t - 5), 



N ■ 



The solution of (llQaD is 



m 



(71 +73|«P) . 



(18a) 
(18b) 

(18c) 



^o(To,Ti) = — + (a(Ti)e^'^^^e^'"°^° +C.C.) , (11 
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The steady state amplitude can be found by setting 
d = 0, nsiow = and taking a square of the evolution 
equation (p!5|) . resulting in 



one obtains (omitting the noise) 



q\l+p^)\af + 2qij,p-Aco)\a\^ 
+ (7? + Ac.2)|a|2 



= 0. (19) 



This cubic equation of |ap can have either one, two, or 
three different real roots, depending on the values of the 
detuning parameter Acj and the excitation amplitude /q. 
When 73 is sufficiently small, i.e., p ^ 0, the solutions of 
(p!9|) behave very much like the ordinary Duffing equation 
solutions, to which ([7|) reduces if 631 = and 632 = (see 
Fig. SI). 




FIG. 4. (Color online) Steady state solutions under different 
excitation amplitudes /o. In case /o < fc (where fc is some 
critical excitation force, dependent on the system parameters, 
see text), only one real solution exists and no bistability is pos- 
sible. In case /o = /c, the system is on the edge of bistability 
and one point exists, where |a|^ vs. to has an infinite slope. 
In case /o > /c, the system is in bistable regime having three 
real solutions over some range of frequencies. Two of these 
solutions are stable. The dashed line denotes the unstable 
solution. 

The solution of (p!5|) can be also presented in polar 
form [41] 



(20) 



where A and (j) are real, and A is assumed to be posi- 
tive. Separating the real and imaginary parts of (p!5|) . 



/o 



A^-fiA^qpA^ = -^sin0, 

ZUJo 

A(j)^ AcoA- qA^ = -— cos(/). 

2cjo 



(21a) 
(21b) 



Steady state solutions are defined by A = 0, ^ = 0, 
which results in ([19]). 

The maximal amplitude |amp can be found from ([19]) 
by requiring 



dAoj 



0, 



where AcOm is the corresponding excitation frequency de- 
tuning. This results in 



AuJn 



3^3 

2uo' 



(22) 



Interestingly enough, the phase (j) of the maximal re- 
sponse is always equal — 7r/2, i.e., the maximal response 
is exactly out of phase with the excitation regardless the 
magnitude of the excitation, a feature well known for the 
linear case. This general feature can be explained as fol- 
lows. For an arbitrary response amplitude A, there exist 
either two or no steady state (j) solutions of ([2T]) . If two 
solutions and ^2 exist, they must obey ^2 = tt — 0i, 



as seen from (|21ap . It follows from (|21bp that these two 
solutions correspond to two different values of Aco. How- 
ever, at the point of maximum response the two solutions 
coincide, resulting in 0i = ^2 = — 7r/2, i.e.. 



(23) 



The system's behavior qualitatively changes when pa- 
rameters such as the excitation amplitude and the fre- 
quency detuning are varied, as seen in Fig. [H The pa- 
rameter values at which these qualitative changes occur 
are called bifurcation (jump) points [46]. 

A jump in amplitude is characterized by the following 
condition: 



or, alternatively. 



dAuj 



dAuj 



±00, 



0. 



Applying this condition to (fT9)) yields 
M\l+p')\aj\'' + ^q{lW-^^j)\aj?- 



(24) 

where Auoj and aj denote the frequency detuning and the 
amplitude at the jump point, respectively. 

When the system is on the edge of bistability, the two 
jump points coincide and ([24|) has a single real solution at 
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the point of critical frequency and critical amplitude 
ttc- The driving force at the critical point is denoted in 
Fig. m as fc- This point is defined by two conditions 



0, 



= 0. 



d{\a\^y 

By applying these conditions one finds 
3^. 



AcJc = y(l+p')|ac|' 



where ac is the corresponding critical amplitude. Substi- 
tuting the last result back into ([24|) . one finds [49.] 



AiJc 71 



271 V3p ± 1 
1 - 3p2 ' 
4p±y3(l+p2) 



1 -3p2 



71 



(26a) 
(26b) 
(26c) 



where the upper sign should be used if 0^3 > 0, and the 
lower sign otherwise. In general, 73 is always positive, 
but as can be either positive or negative. Therefore, q 
and p are negative if a3 < (soft spring) , and positive if 
as > (hard spring). 

It follows from (|26ap that the condition for the critical 
point to exist is 



bl< 



Without loss of generality, we will focus on the case of 
"hard" spring, i.e., 0^3 > 0, g > 0, p > 0, as this is the 
case encountered in our experiments. 



C. Behavior near bifurcation points 

When the system approaches the bifurcation points, it 
exhibits some interesting features not existent elsewhere 
in the parametric phase space. In order to investigate the 
system's behavior in the vicinity of the jump points, it 
is useful to rewrite the slow envelope evolution equation 
(p!5|) as a two dimensional fiow 



(27a) 
(27b) 



where we have defined x(t) = Re{a}, y{t) = Im{a} (i.e., 

a{t) = x{t) -\- jy{t)), and 



g{x, y) = -Aujx - jiy + q{x^ + y^){x - py) 



y^){y^px), (28a) 

A. 

2cjo* 
(28b) 



The real- valued noise processes nx{t) and ny{t) have the 
following statistical properties: 



K(t)) = (n,(t)) =0, 
K(t)n,(t)) =0, 

{ncit)nc{s)) = {nsit)nsis)) 



(29a) 
(29b) 



^S{t-s). (29c) 



At the fixed points, the following holds: f{x^y) = 
g{x,y) = 0. A typical phase space fiow of the oscilla- 
tor in bistable regime is shown in Fig. [H 
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FIG. 5. (Color online) Slow envelope phase plane trajectories 
of a nonlinear oscillator in bistable regime. Three real solu- 
tions of (p!9|) correspond to three fixed points of the flow. 5*1 
and S2 are stable spiral nodes, whereas U is the saddle-node, 
from which two manifolds emerge [46]. The green dotted line 
is the stable manifold ("separatrix"), which separates differ- 
ent basins of attraction, belonging to different stable nodes 
Si and S2. The magenta thick line is the unstable manifold. 
Two typical phase plane trajectories are shown by arrowed 
thin blue lines. 

For small displacements near an arbitrary fixed point 
ao = (xo,^o), namely, x = xq ^ Ax and y = yo ^ Ay, 
where Ax <C xq and Ay <C yoj the above nonlinear fiow 
map can be approximated by its linearized counterpart 



Ax 
Ay 



HZ 



where 



M 



fx fy 
9x 9v 



(30) 



(31) 



and the excitation frequency detuning Auj as well as the 
external excitation amplitude /o are considered constant. 
The subscripts in the matrix elements denote partial 
derivatives evaluated at (xo,^o), for example, 

dx x=xo, 

y^yo 
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The matrix M is, therefore, the Jacobian matrix of the 
system ([27|) evaluated at the point (xq, ^o)- It is straight- 
forward to show that 



fx = -7i - QVi^l + yl) - '^QXoiVo +P^o), 
fy = Au;- q{xl + y'^) - 2qyo{yo ^ pxo), 
gx = -Aij + q{xl + yl) + 2qxo{xo - pyo), 
9y = -71 - QPi^l + Vo) + 2qyo{xo - pyo). 

Two important relations follow immediately, 

f = 9yX- gxV, 



(32a) 
(32b) 
(32c) 
(32d) 

(33a) 
(33b) 

The linearized system (|3Q|) retains the general quali- 
tative structure of the flow near the flxed points [40|, in 
particular both eigenvalues of the matrix M are negative 
at the stable nodes, denoted as Si and ^2 in Fig. \S\ At 
the saddle node [/, which is not stable, one eigenvalue of 
M is positive, whereas the other is negative. 

The discussed Dufling like systems exhibit saddle point 
bifurcations. At the bifurcation, one of the stable nodes 
and the saddle node coincide, resulting in a zero eigen- 
value in M. The bifurcation ("jump") point condition is, 
therefore, det M = 0, which gives the same result as in 
([24|) . The case of well separated stable node and saddle 
node is shown in Fig. [6l and the case of almost coincid- 
ing stable and not stable flxed points is shown in Fig. [3 
where the oscillator is on the verse of bifurcation. 
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FIG. 6. (Color online) The phase plane geometry when the 
saddle-node (U) and the stable node (6*2) are well separated. 
The green dotted line is the stable manifold ("separatrix") 
and the magenta thick line is the unstable manifold. A typi- 
cal phase plane trajectory is shown by the arrowed thin blue 
line. The absolute value of slow envelope's rate of change a 
is represented by the background color. The darkest parts 
correspond to the slowest motion in the phase space. At both 
fixed points U and 6*2 the value of d is zero. 
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FIG. 7. (Color online) The phase plane geometry when the 
saddle-node (U) and the stable node {S2) are close one to 
another. The green dotted line is the stable manifold ("sepa- 
ratrix"), and the magenta thick line is the unstable manifold. 
Phase plane trajectories are shown by the thin blue lines. The 
absolute value of d is represented by the background color. 
The darkest parts correspond to the slowest motion in the 
phase space. At both fixed points U and &, the value of d 
is zero. Note that the motion in the phase space becomes es- 
sentially one-dimensional and slows down significantly in the 
vicinity of the stable node 6*2 . 



We note that in general 
TrM = /,+^^ = 



-2(7i + 2gp|aop), 



and the slow eigenvalue near the bifurcation point can be 
estimated as 



detjd d det M 
Xr M ^ dAuj TrM 
V 



■5 



Auj=Auj 
a—aj 

2q\aj 



71 + 2qp\aj 



where (5 is a small frequency detuning from Acjj, i.e. 
Aco = Aujj + ^, 1^1 <C |Acjj|. If the system in bistable 
regime is close to bifurcation then Agd and the evo- 
lution of the system almost comes to stagnation, phe- 
nomena often referred to as critical slowing down [49|. 
The motion in the vicinity of the stable node becomes 
slow and essentially one-dimensional along the unstable 
manifold. We now turn to show this analytically. 

At the bifurcation points the matrix M is singular , 
i.e., det M = 0. Consequently, the raws of the matrix are 
linearly dependent, i.e.. 



^ I fx fy 

rfx rfy 



where r is some real constant. Using the last result, we 
may rewrite (|33ap at the bifurcation point as 



r {fyx - fxV) = 0, 
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where we have used the fact that at any fixed point 
(stable or saddle-node) f{x^y) = g{x,y) = 0. How- 
ever, according to (|33bp . at any fixed point fyX — fxV = 
—fo/2uJo 7^ 0. Therefore, r = at the bifurcation point, 
and the matrix M can be written as 



M = A/ 



1 K 
I 



where 



A/ 
K -- 



:-2(7i + 2gp|a,f ) 
fy _ 7i -^2qp\aj 



fx 2q\aj\^ ■ 



It also follows from (l33aD that 



71 ^QP\aj\" 
Aoo^Aujj gx AcOj — q\aj\' 



^= lim ^ 



(34) 

(35a) 
(35b) 

(36) 



Due to the singularity of matrix M at the bifurcation 
point, a second order Taylor expansion must be used. 
The fiow map ([27|) can be approximated near the bifur- 
cation point by 



Ax = Xf{Ax + KAy) + fsS 



Ay ggS 
1 

+ 2 



d 



dAuo 



Ax 



d_ 

dx 



Ay 



dy 



/ + n^(t), (37a) 



dAcj 



Ax 



dx 



Ay 



dy 



g^ny(t), (37b) 



where all the derivatives denoted by subscripts are eval- 
uated at the jump point a = aj, and 

gs -Xj. 

The above system of differential equations (|37|) can be 
simplified by using the following rotation transformation, 
shown in Fig. [H 



cos a sm a 
— sin a cos a 



Ax 
Ay 



(38) 



where tana = K. In these new coordinates, the system 
(j57)l becomes 

^ = Xf^ + n^5+]^B^H{^,r,)+n^{t), (39a) 
y? = -XfK^ + 0^.5 + ^ E{^, r?) + n^(t), (39b) 



where 



v) = /(^ , V) cos a + g{^, rj) sin a, 
^(^, V) = 9(^1 V) cos a - /(^, T]) sin a, 
= yj cos a — Xj sin a. 



^71 



-Xj COS a — yj sma. 



and D is the differentiation operator 



d_ 



V 



= Ax 



dx 



d_ 
df] 

-Ay 



d 



dAuj 

d , d 
k A — 

dy 



dAuj 



The noise processes n^{t) and n^(t) have the same sta- 
tistical properties ([29|) as nc(t) and ns{t). 
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FIG. 8. The effective one dimensional flow in the vicinity of 
a bifurcation point. The origin of the phase plane coincides 
with the bifurcation point. U is the saddle point, and S2 is 
a stable node. The effective flow between these two points, 
marked by arrows, is almost parallel to the rotated coordinate 
?7, while the rotated coordinate ^ remains essentially constant, 
^ = ^0 + 0(?7^). a IS the angle of coordinate rotation. The 
velocity of the flow is largest at the point M, between the 
saddle point and the stable node. 

The time evolution of the system described by the dif- 
ferential equations (|39|) has two distinct time scales. Mo- 
tion along the coordinate ^ is "fast", and settling time 
is of order |A/|~^. The time development along the co- 
ordinate 7^, however, is much slower, as will be shown 
below. 

On a time scale much longer than |A/|~^, the coordi- 
nate ^ can be regarded as not explicitly dependent on 
time. The momentary value of ^ can be approximated as 



1 



1 d^H 

2 Ot]'' 



-V 



(40) 



where we have neglected all terms proportional to S'^ and 
St]. 

The motion along the coordinate r] is governed by a 
slow evolution equation (|39bp , combining which with (|4Q|) 
results in 



77 = 7)0 Bt]'^ -^rirjit), 



where 



^0 



(41) 



(42a) 
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B = — - — \xj(l + 2 sin^ a -\- p sin 2a) 
cos a ^ 

- yj + 2 cos^ a) + sin 2a)] . (42t 

Note that the noise process n^{t) does not play a si^ 
nificant role in the dynamics of the system, because th 
system is strongly confined in ^ direction. Such nois 
squeezing is a general feature of systems nearing saddk 
point bifurcation [25l, [H, [H, (SS^ . 

Two qualitatively different cases of (|4T]) should be rec 
ognized. The first case is of a system in a bistable regim 
with a stable (quasi stable, as we will see below) and no 
stable (saddle) fixed points close enough to a bifurcatio 
point. In this case, the one dimensional motion is equi\ 
alent to a motion of a massless particle in a confinin 
cubic "potential" 

U{v) = -V (vo + isr?') , (43) 
as shown in Fisr. [9l Fisrure [10] deoicts the location of 
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FIG. 9. (Color online) Effective one dimensional potential 
U{r]) oc -ry (770 + |Bry2) i^. 

the fixed points and the bifurcation point on a frequency 
response curve in this case. Figure [TT] shows a comparison 
between the exact simulation of the system's motion near 
the bifurcation point and the analytical result (|4T]) . 

The quasi one-dimensional system described above is 
obviously not stable j43,[8l|. The rate of escape from the 
vicinity of the quasi stable fixed point is 0,l85| 

rtheTm[0) ^ — e D ^ 




Aujj 



FIG. 10. (Color online) The location of stable nodes Si and 
5*2, and a saddle node [/ in a bistable regime close to a bi- 
furcation point. 6^ which is negative in this case, is the fre- 
quency difference between the excitation frequency and the 
jump point frequency (cjo + Acjo) + Acjj. The scales of the 
axes are arbitrary. 
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FIG. 11. (Color online) Velocity along the slow coordinate 
r\ for different values of detuning 8. Aujj is the jump point 
(bifurcation) detuning, p = 0.3/ VS in all cases, U, U' and U" 
are the saddle node positions for different values of 5. Simi- 
larly, 5*2, S2 and S2 are the stable node positions for different 
values of 5. Exact values of drj/dt are shown by solid lines. 
The dashed lines are the results of analytical approximation 
(|4T]) . The scales of the axes are arbitrary. 
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where 



D 



N 



(71 +73|%f 



stable node 



Characteristic time of thermal escape rthe 
shown to be [42i] 



'H^hermC^) — 



Hher 



Toe 



A 

fcRT 



where 



To 



TT 



64 



(X 



i-sy 



3 71 +73|«j| 



can be 
(44a) 

(44b) 
(44c) 



This is a mean time in which the system escapes from 
the stable node near bifurcation point to the other stable 
solution of (p!5]) due to thermal noise nrj{t), and the 3/2 
power law is correct as long as A ksT [48]. 

The second case describes a system which has under- 
gone saddle bifurcation, i.e., an annihilation of the stable 
and non stable points has occurred. The phase plane 
motion close to the bifurcation point is still one dimen- 
sional, however, 7)0 changes its sign. Therefore, the mo- 
tion is not confined any more, but is still very slow in 
the vicinity of the bifurcation point, because fjo (x S, as 
follows from Eq. (|42ap . The system starts converging 
to the single remaining stable fixed point, but is signif- 
icantly slowed down, and lingers in the vicinity of the 
bifurcation point due to the saddle node "ghost". As the 
system spends most of its time of convergence near the 
saddle node "ghost", this slow time of convergence r^d 
can be roughly estimated as [46|] 



f 

Jo 



TT 



(45) 



Note that Tgd cx due to (|42ap . 

D. Extraction of parameters from experimental 
data 

The analytical results presented above allow us to use 
data acquired in relatively simple experiments in order 
to estimate several important dynamic parameters of the 
micromechanical beam. We note that data acquisition 
using e-beam or optical beam interaction with vibrating 
elastic element does not readily enable extraction of dis- 
placement values. In contrast, the frequencies of impor- 
tant dynamical features, including maximum and jump 



points, can be measured with high accuracy using stan- 
dard laboratory equipment, such as network analyzers 
and lockin amplifiers. Therefore, it is desirable to be able 
to extract as much data as possible from the frequency 
measurements. 

If the system can be brought to the verge of bistable 
regime, i.e., /o = /c, the nonlinear damping parameter p 
can be readily determined using Eq. (|26cp . The same co- 
efficient can also be extracted from the measurements of 
the oscillator's frequency response in the bistable regime. 
In general, the sum of the three solutions for |ap at any 
given frequency can be found from Eq. (p!9|) . This is em- 
ployed for the jump point at loq + Acoq + Acoj seen in Fig. 
m Using Eq. ([22|) to calibrate the measured response at 
this jump point one has 



(2/li + /l2)|am|' = 



2q (71P - Aujj 



or 



{2hi + h2) Aujmil^p'^) + 2 (71P - Au 



0, (46) 



where hi and /12 are defined in Fig. [H Due to the fre- 
quency proximity between the maximum point and the 
jump point at = cjo+Acjo+Acj^, the inaccuracy of such 
a calibration is small. Moreover, as long as excitation am- 
plitude is high enough, /12 is much smaller than hi and 
even considerable inaccuracy in /12 estimation will not 
have any significant impact. This equation can be used 
to estimate p for different excitation amplitudes at which 
the micromechanical oscillator exhibits bistable behav- 
ior, i.e., /o > fc- It is especially useful if the system is 
strongly nonlinear and cannot be measured near its crit- 
ical point due to high noise ffoor or low sensitivity of the 
displacement detectors used. 

Another method for estimating the value of p requires 
measurement of free ring down transient of the microme- 
chanical oscillator and can be employed also at low ex- 
citations, when the system does not exhibit bistable be- 
havior, i.e., /o < fc- The polar form of the evolution 
equation ([2T]) is especially well suited for the analysis of 
the system's behavior in time domain. Starting from 
Eq. (|21ap and applying the free ring down condition 
/o = 0, one finds 



Ale- 



1 



qp 
71 



-271^ 



(47) 



where Aq is the amplitude at t = 0. In particular, con- 
sider a case in which the system is excited at its maxi- 
mal response frequency detuning Acj^, i.e., Aq = |amp. 
Then, after turning the excitation off, the amplitude dur- 
ing the free ring down process described by Eq. (|47|) can 
be written as 



AHt) 



-27it 



-p 



(l-e-27it) 



(48) 



The ring down amplitude measured in time domain can 
be fitted to the last result. 
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In addition to nonlinear damping parameter p, most 
parameters defined above can be easily estimated from 
frequency measurements near the jump point shown in 
Fig. [4] if the following conditions are satisfied. The first 
condition is 



«1, 



(49a) 



which can be satisfied by exciting the micromechanical 
beam oscillator in the bistable regime strongly enough, 
i.e., /o ^ fc- The immediate consequence of the first 
condition is 



«1, 



As shown in Sec. IIIIB[ Eq. ([23|) . at the maximum re- 
sponse point Au = Aumj the following holds: = 
—j\am\- Therefore, in view of our assumptions described 
above, we may write 



Consequently, 



(49b) and 



^0 



(51a) 



i.e., hi ^ 1, diS described above. 

Using Eq. (|49bp . it follows from Eq. lucil (^|u.^ | 
AcUjri- From the last result and from Eqs. (|35|) . (|36|) . and 
(|38|) . the following approximations follow immediately: 



that ^1%'^ 



K = tan a 
1 . 



71 

AcOn 



2p, 



A/ 



K ' 

-2(7i+2pAcj^), 
7i + pAum 



Aujj 



AuJn 



(50a) 

(50b) 
(50c) 
(50d) 



B 



AUn 



\aj\{l + K^y. 



-I 



p (3 + if2) 



2K 



(51b) 



which follows from Eq. (|12|l . The time Tgd, which de- 
scribes the slowing down near the saddle-node "ghost" 
described above Eq. can be expressed as 



where 



Y = 



Tsd((5) 



271' 




r 



(52) 



(53) 



Finally, we turn to estimate the value of the thermal 
escape time rtherm given by Eq. (|44|) . Using the same 
assumptions as above, we find 



To 



A/7 



64 



(54a) 



3 7l +73|«ri 



(54b) 



Unlike in the previous approximations, one has to know 
at least the order of magnitude of the response amplitude 
in the vicinity of the jump point (in addition to effective 
noise temperature T and effective mass m) in order to ap- 
proximate Ttherm appropriately. The same is also true for 



estimation attempt of the physical nonlinear constants 

2ujQAuJrn 



as 



73 2p- 



3|amp 

AuJn,^ 



12- 



(55a) 
(55b) 



For more accurate estimation, one of several existing 
kinds of fitting procedures can be utilized tSS, £1]. How- 
ever, the order of magnitude estimations often fully sat- 
isfy the practical requirements. 
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E. Experimental considerations 

The above discussion of parameters' evaluation using 
experimental data, especially in frequency domain, em- 
phasizes the importance of accurate frequency measure- 
ments. However, the slowing down of the oscillator's re- 
sponse near the bifurcation points poses strict limitations 
on the rates of excitation frequency or amplitude sweeps 
used in such measurements [86|. This is to say that spe- 
cial care must be taken by the experimentalist choosing 
a correct sweep rate for the measurement in order to ob- 
tain the smallest error possible. Fortunately, this error 
can easily be estimated based on our previous analysis. 

Let Tsweep represent the frequency sweep rate in the fre- 
quency response measurement. For example, using net- 
work analyzer in part of our experiments, we define 



' sweep 



27r 



frequency span ( Hz) 
sweep time (sec) 



(56) 



In order to estimate the inaccuracy, (5err, in the measured 
value of the bifurcation point detuning, Auj, which re- 
sults from nonzero frequency sweep rate, the following 
expression may be used: 



\Se 



I sweep 



'^sd (^err)? 



whose solution is 



(-Yr V 

^2 sweep j 



(57) 



Note that this error is a systematic one - the measured 
jump point will always be shifted in the direction of the 
frequency sweep. Obviously, the first step towards accu- 
rate measuring of Acjj is to ensure that the established 
value of the bifurcation point detuning does not change 
when the sweep rate is further reduced. 

Another possible source of uncertainty in frequency 
measurements near the bifurcation point is the thermal 
escape process. The error introduced by this process 
tends to shift the measured jump point detuning in the 
direction opposite to the direction of the frequency sweep. 
Moreover, unlike the error arising from slowing down pro- 
cess, this inaccuracy cannot be totally eliminated by re- 
ducing the sweep rate. However, as will be shown in 
Sec. IIVB[ in our case this error is negligible. 



and damping backbone curve depicted in Figs. [T3l and [Ml 
respectively, for a 125 pm long beam with fundamental 
mode resonant frequency 885.53 kHz and Vdc = 20 V. 
We derive the value of 71 = uJo/2Q from the linear re- 
sponse at low excitation amplitude and find Q = 7200 
for 200 pm beam and Q = 13600 for 125 pm beam. 
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Frequency (kHz) 



123.39 123.45 



FIG. 12. (Color online) Measured response amplitude vs. 
excitation frequency shown for both upward and downward 
frequency sweeps with Vdc = 20 V and with varying peak- 
to-peak excitation amplitude Vac of a 200 pm long beam with 
fundamental mode occurring at 123.2 kHz. The excitation 
amplitude is shown on the graphs. The oscillator exhibits 
bistable behavior at all excitation amplitudes except for the 
lowest one. The vertical axis is in arbitrary units. 



As shown in Sec. lIIIDl the value of p can be estimated 
for different excitation amplitudes using Eqs. (|46|) and 
(|48|) . Typical results of applying these methods to ex- 
perimental data from a micromechanical beam oscillator 
can be seen in Fig. [151 Using these procedures we find 
p ^ 0.292 for the 200 pm long beam and p ^ 0.109 for 
the 125 pm long beam. We also estimate p from the crit- 
ical point detuning using Eq. (|26cp . and obtain similar 
values. 



IV. RESULTS 

A. Nonlinear damping 

A typical measured response of the fundamental mode 
of a 200 pm long beam occurring at the resonance fre- 
quency of 123.2 kHz measured with Vdc = 20 V and 
varying excitation amplitude is seen in Fig. [121 The lin- 
ear regime is shown in the frequency response diagram 



B. Parameter evaluation 

In order to illustrate the procedures derived in 
Sec. IIIID[ we evaluate the main parameters of slow en- 
velope dynamics of the 125 pm long beam in a particular 
case in which Vdc = 15 V and the excitation voltage 
amplitude is 140 mV. The quality factor of the beam, as 
measured in the linear regime, is Q = 13600. 
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FIG. 13. (Color online) Measured response vs. frequency in a 
linear regime of the 125 pm long beam with fundamental mode 
occurring at 885.53 kHz and Vdc = 15 V. The linear regime is 
defined as a regime in which the frequency response function 
is symmetric around the resonance frequency. In the main 
panel, the measured responses with three different excitation 
amplitudes are shown. Blue circles correspond to ^; = 10 mV, 
green rectangles correspond to v = 20 mV, and red triangles 
correspond to ^; = 30 mV. Solid black lines show the fitted 
Lorentzian shapes. Vertical scale is in arbitrary units. 
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FIG. 14. (Color online) Measured response amplitude \a\ vs. 
inverse quality factor 1/Q in a linear regime [3, [g^I. Large 
black diamonds correspond to the frequency responses de- 
picted in Fig. [131 The measured averaged quality factor is 
Q = 13600 ± 4%. Other experimental parameters are simi- 
lar to those described in Fig. [13] caption. Vertical scale is in 
arbitrary units. 



The results that can be derived from frequency 
measurements only, i.e., the results corresponding to 
Eqs. (|26cp . (f5Ql) and ([53l), are summarized in Table [D 

For this measurement we employ a network analyzer 
with frequency span of 500 Hz, sweep time of 13.6 sec, and 
bandwidth of 18 Hz. Therefore, the sweep rate defined in 
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FIG. 15. (Color online) Experimental results for p = 
73a;o/3a3 vs. excitation amplitude. The excitation ampli- 
tude on the horizontal axis is normalized by the respective 
critical excitation amplitude fc. a. 200 pm long beam with 
fundamental mode occurring at 123.20 kHz and Q = 7200. 
The values of p extracted from frequency domain jump point 
measurements (see Eq. (|46|) ) are represented by blue circles. 
Red dashed line represents the value oi p — 0.505/ VS = 0.292 
evaluated using the critical point frequency detuning AcOc 
(see Eq. ()26cp ). The critical excitation voltage is 50 mV, and 
Vdc = 20 V. b. 120 pm long beam with fundamental mode 
occurring at 885.53 kHz and Q = 13600. The values of p 
extracted from time domain ring down measurements accord- 
ing to Eq. (|48|) and frequency domain jump point measure- 
ments (see Eq. ()46p ) are represented by green squares and blue 
circles respectively. Red dashed line represents the value of 
p = 0.189/a/3 = 0.109 evaluated using the critical point fre- 
quency detuning AcOc that is given by Eq. ()26cp . The critical 
excitation voltage is 105 mV, and Vdc = 15 V. 



Parameter 


Value 


Units 


a;o/27r 


885534 


Hz 


71 


204 


sec""*^ 


AoOm/^TT 


76 


Hz 


AiOj /27r 


81 


Hz 


P 


0.109 




K 


0.646 




a 


0.573 


rad 


A/ 


-616.5 


sec""*^ 


y± 

Xj 


8.16 




Y 


0.158 


1 

sec 2 



TABLE I. Parameters of the slow envelope dynamics of a 
125 pm long beam. Applied DC voltage is 15 V and excitation 
voltage amplitude is 140 mV. The critical excitation voltage 
is 105 mV. Quality factor is Q = 13600. 



Eq. ([561) is 



' sweep 



o 500 Hz ^ _2 

ZTT =231radsec 

13.6 sec 
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The inaccuracy in jump point detuning estimation due 
to slowing down process (see Eq. ([57|) ) is 

|^^2Hz. (58) 

ZTT 

We now turn to estimate the order of magnitude of 
other parameters, including the nonlinear elastic con- 
stant 0^3 and nonlinear damping constant 73. Based on 
the observations of the vibrating micromechanical beam 
by the means of SEM continuous scanning mode, we esti- 
mate the amplitude of mechanical vibration to be around 
100 nm. The mass of a golden beam of the dimensions 
given in Sec. [Til is approximately 7 x 10~^^kg. These es- 
timations allow us to assess the order of magnitude of 
several additional parameters shown in Table [Til which 
is based on Eqs. ([54|) and ([55|) . 



Parameter 


Value at 


Units 


6 = 


— 6err — — 27r X 2 Hz 




as 


2 X 10^^ 


— 2 —2 

m sec 


73 


1 X 10^^ 




T 


300 


°K 


AU/ksT 


6 X 10^ 




To 


0.13 


sec 



TABLE II. Order of magnitude estimation of parameters of 
a 125 pm long beam's slow envelope dynamics. The distance 
from the excitation frequency to the jump frequency is taken 
to be equal to ^err (see Eq. ([58]) ). Applied DC voltage is 15 V, 
the excitation voltage amplitude is 140 mV, and the estimated 
amplitude of vibration is 100 nm. The critical excitation volt- 
age is 105 mV. Quality factor is Q = 13600. 

We estimate below the thermal escape time for 5 = 
— Jerr (scc Eq. ([58]) ). Howcvcr, the value of the expo- 
nent, AU/ksT - 6 X 10^ at T = 300 K, makes the ther- 
mal escape time at this detuning value extremely large. 
Therefore, in our experiments, the thermal escape pro- 
cess does not contribute significantly to the total inac- 
curacy in frequency measurements near the bifurcation 
point, at least for effective noise temperatures lower than 
10^ K, at which the assumption A ksT is no longer 
valid. 

Finally, it is interesting to compare the nonlinear dis- 
sipation term 73|ap and the linear dissipation term 71 in 
the evolution equation ([T5|) . It follows from the above as- 
sumptions and the values in Table [Til that for our chosen 
example 



71 

C. Validity of the multiple scales approximation 

In order to verify the correctness of our approximated 
solution achieved by multiple scales method, we compare 



the results of direct integration of the full motion equa- 
tion ([9]) with the steady state solution of the evolution 
equation ([T9|) . We use the results from Tables [Tj and [Til 
for c^o, <^3, 71, and 73. We also estimate the effective mass 
m to be 0.7 X 10~^^ kg, the effective capacitance to be of 
order of Co ^ 1.5 x 10"^^ F, the DC voltage Vbc = 15 V, 
the AC voltage v = 200 mV, and take the distance d 
to be the actual distance between the electrode and the 
mechanical beam, i.e., d = 5 pm. The resulting exci- 
tation force amplitude is /o = 600 Nm~^, the constant 
force is F = 45000 Nm"^ (see Eq. and the constant 
resonance frequency shift is Acoq = — 27r x 257 Hz (see 
Eq. ([T4|)). 

In Fig. [TBI the exact numerical integration of Eq. ([9]) 
is compared with the solution of the approximated fre- 
quency response equation ([T9|) . A very good correspon- 
dence between the two solutions is achieved, which vali- 
dates the approximations applied in Sec. IIIIBi 
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FIG. 16. (Color online) Comparison of numerically calculated 
steady state response amplitude of the full equation of motion 
(|9]) (red circles) with the steady state solution of the evolution 
equation ([19)) (solid line). 



V. DISCUSSION 
A. Analysis of results 

It follows from our experimental results that the non- 
linear damping constant p can be estimated with a high 
degree of confidence by measuring the micromechanical 
oscillator bistable response in the frequency domain. The 
values of p that we find, 0.1 < p < 0.3, obviously are 
not negligible. Referring to Eqs. (|26ap and ([24|) . we see 
that the considered micromechanical oscillators exhibit a 
damping nonlinearity that has a measurable impact on 
both the amplitude and frequency offset of the critical 
point, as well as on jump points in the bistable region. 
On the other hand, these values are significantly smaller 
then the critical value p = 1/a/3 ^ 0.577, which would 
prevent the system from exhibiting bistable behavior. 

Two methods of estimating the value of p from fre- 
quency domain measurements were used. The first is 
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based on a single measurement of the critical point and 
provides a simple means for estimating the value of p by 
experimentally measuring the linear quality factor Q at 
low excitation amplitude and the critical frequency shift 
AuJc only (see Eq. (|26cp ). The second can be used for any 
excitation amplitude that drives the system into bistable 
regime, but requires a comparison of different response 
amplitudes (see Eq. (|46|) ). Both these methods yield sim- 
ilar results, however, the second one, although being less 
accurate, allows the experimentalist to estimate when the 
limit of hard excitation [41] is approached and the first 
order multiple scales analysis used in this study becomes 
inadequate. In this limit of strong excitation, the ex- 
tracted values of p start to diverge significantly from the 
results obtained at low excitation amplitudes. Our re- 
sults, especially Fig. [151 ^nd the analysis of the validity of 
our approximations, which was carried out in Sec. IIV 
suggest that the analysis method employed by us is ade- 
quate for a wide range of excitation amplitudes. 

The third method described above allows one to es- 
timate the value of p from time domain measurements 
of the free ring down of the micromechanical beam os- 
cillator based on Eq. (|48|) . Although fitting results of 
time domain measurements to a theoretical curve intro- 
duces large inaccuracy, this method is invaluable in cases 
where the bistable regime cannot be achieved, e.g., due 
to prohibitively large amplitudes involved and the risk of 
pull-in. 

By using the approximations developed in Sec. IIIID[ 
we were able to estimate different parameters describing 
the slow envelope dynamics of our oscillators, summa- 
rized in Table H The most important and, as far as we 
know, novel result is the direct estimation of the slowing 
down time r^d that is given by Eq. ([52|) . which governs 
the system's dynamics in the vicinity of bifurcation point. 
In turn, this result is used to quantitatively evaluate the 
error introduced to the frequency measurements by the 
slowing down process, ^err that is given by Eq. ([57|) . which 
in the example studied is 2 Hz. It can be seen that even 
slow sweeping rate (as compared to quasistatic rate in the 
linear case, which is of order of one resonant width per 
ring down time) can introduce a significant inaccuracy in 
the measured response of a micromechanical beam oscil- 
lator near bifurcation points. In our case, the inaccuracy 
in Aujj is about 3%, but the inaccuracy in Auoj — AuOm 
is probably much larger. 

The nonlinear damping constant p plays an important 
role in all the dynamical parameters. In the value of K 
that isgiven by Eq. (|5Qap in our example, p-dependent 
term constitutes about 30% of the value. The same is 
true for other parameters as well. 

Also, we make order of magnitude estimations of ther- 
mal escape time rthermai (see Eq. ([54|) ). 0^3, and 73 (see 
Eq. ([55]) ). which are summarized in Table HIl These ap- 
proximations can be used in order to construct an ac- 
curate model of the effective one-dimensional movement 
of the system in the vicinity of a bifurcation point, es- 
pecially if accurate enough estimations of the oscillator's 



amplitude and effective mass can be made. 

In our case, only the order of magnitude of the parame- 
ters can be estimated. However, we were able to estimate 
the thermal escape time, and found the thermal escape 
process to be a non negligible source of inaccuracy in the 
frequency measurements only at very high effective noise 
temperatures of order 10^ K. This result can be com- 
pared to a result from our previous work [25]. In that 
work, a micromechanical beam oscillator similar to the 
ones used here was excited at a frequency between the bi- 
furcation points. The intensity of voltage noise needed to 
cause transitions between these stable states was found to 
be ^ 500 mV, with noise bandwidth of 10 MHz. The re- 
sulting voltage noise density is 158 /iV/V^Tlz, which cor- 
responds to an effective noise temperature ~ 10^^ K. In 
the case of thermal escape described here, the two stable 
states are highly asymmetrical. The effective noise tem- 
perature of 10^ K, which invalidates the estimations of 
very slow thermal escape rate in Sec. lIVBt corresponds 
to voltage noise density of 0.5 /iV/VHz, giving the total 
voltage noise intensity of 1.6 mV. 

Finally, we can also estimate the relative contribution 
of the nonlinear damping term 73|ap in the evolution 
equation (p!5]) . and find it to be a non negligible one- 
tenth of the linear term 71 (see Eq. ([59]) ) at the estimated 
amplitude of \a\ = 100 nm. 

B. Geometric nonlinearities as a source of 
nonlinear damping 

The nature of nonlinear damping is not discussed in 
this work. However, nonlinear damping can be, in part, 
closely related to material behavior with a linear dissi- 
pation law that operates within a geometrically nonlin- 
ear regime. Here, we investigate one possible mechanism, 
originating from a Voigt-Kelvin type of dissipation model 
which describes internal viscoelastic damping in the form 
of a parallel spring and dashpot. 

Before we proceed to build the model, one technical re- 
mark is in order. The notations in this section follow the 
standard ones used in continuum mechanics, and some 
parameters used above are redefined below. However, 
the end results are brought back to the form of ([8]) . 

Following Leamy and Gottlieb [s^, [H| , we consider a 
planar weakly nonlinear pretensioned, viscoelastic string 
augmented by linear Euler-Bernoulli bending, which in- 
corporates a Voigt-Kelvin constitutive relationship where 
the stress is a linear function of the strain and strain rate 
|6i,[8i|: 

a = Ee^ Det, 

where a is the stress, e is the strain, E is the material 
Young modulus, I) is a viscoelastic damping parameter, 
and subscripts denote differentiation with respect to the 
corresponding variable. The equations of motion of the 



17 



beam-string are 



pAutt 



Nus + EA ( Us + -wl 



DA {Uts + WsWts) 



0, (60a) 



Note that defines the ratio between the longitudinal 
and transverse wave speeds [4l|,[8§|. The rescaled parallel 
plate approximation is thus: 



(7 — wy 



where 



pAwtt 



Nws + EAws I Us + -^^ 



-DAWs {uts + WsWts) - {EIWsss + DlWtsss) 



(60b) 



where N is the pretension, p is the material density, 
s is the material coordinate along the beam, A and / 
are the elastic element cross-sectional area and moment 
of inertia, respectively. Also, u{s^t) and w{s^t) are the 
respective longitudinal and transverse components of an 
elastic field. The generalized transverse force component 
Qu, is due to external electrodynamic actuation. Note 
that for a parallel plate approximation. 



B 



{d-wY 



where Fdc, ^, d and uo are as those defined in Eq. (|8]), 
and 5 is a proportionality coefficient dependent on the 
exact geometry of the mechanical oscillator. 

We rescale the elastic field components u and and 
the material coordin ate s by th e beam length L, and time 
by the pretension ^J pAL?' jN to yield a coupled set of 
dimensionless partial differential equations for the beam- 
string: 



Us^ ^\Us^ -^s ) + (5 [Urs + y^sy^Ts) 



= 0, 
(61a) 



Ws + t^Ws ( Us + -1^5 ) + ^Ws {Uts + WsWrs) 



- {aWsss + pSWrsss) 



where u = u/L, w = w/L, s = s/L and 



Qw. (61b) 



V 



LN ' 



UJs 



We note that as the first longitudinal natural frequency 
is much higher than the first transverse natural frequency 
^ 1), the longitudinal inertia and damping terms 
in Eq. (|61ap can be neglected to yield a simple spa- 
tial relationship between the transverse and longitudi- 
nal derivatives. Incorporating fixed boundary conditions 
(ii(0,r) = u{l^r) = 0) enables integration of the result- 
ing relationship to yield: 



where 



1 2 

2 ' 



•ci(r), 



wzds. 



Thus, the resulting weakly nonlinear beam-string initial 
boundary value problem consists of an integro-differential 
equation for the transverse mode: 



Wr 



where 



Wss (1 + Pci{t) + Scir{r)) 

+ aWssss + pSWrs 



Qw, (63) 
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WsWrsds. 



In order to facilitate comparison of the continuum model 
with the lumped mass model in Eq. (|8]), we consider a 
localized electrodynamic force Qw = Qw{s = 1/2, r). 

We reduce the integro-differential field equation in (|63|) 
and its fixed boundary conditions to a modal dynamical 
system via an assumed single mode Galerkin assumption, 
w{s,t) = qi{r)(j)i{s), using a harmonic string mode = 
^/2sm{7^s): 



N 



' pAL^ ' 

Other dimensionless parameters include the effects of 
weak bending a < 1, a strong nonlinear pretension > 1, 
a small slenderness ratio p < 1 (because r/L <C 1, where 
r = \p~j~A is the beam-string radius of gyration [89]), 
and finite viscoelastic damping ^: 



h\\i- 



SqQr 

[1 



+ h {oiQ + pSqr) 



(7 - hqY 



(64) 
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where q = qi and the integral coefficients are: 



Ii= [ (t)\ds = 1 
Jo 

f 

Jo 



h 
h 
h 
h 



(pi4>usds = -n , 



i>i4'issssds = n 



It is convenient to rescale the maximal response 
\w{l/2,r)\ = q4>, where_^ = 0i(l/2) = V^, by the di- 
mensionless gap z = qcl^/jj and to rescale time by the 
unperturbed {r] = 0) natural frequency t' = cDir, where 
uji = \J al4 — I2 = ttVI + Q^TT^. The resulting dynamical 
system is: 



where 



^2U^^[l+eCOs(m')F 



(65) 



o_ 1 W 

1 I4jii5 

Q Cbi 

Note that the ratio between nonlinear and linear damping 
in Eq. (|65|) consists of only the beam-string geometric 
properties jsO]. For example, a typical ratio is SQ = 
Qd? /h? ^ 65 for a beam-string with a prismatic cross- 
section, where h = 1.5 pm is the dimension of the beam- 
string in the transverse direction and d = 5 pm is the 
resonator gap. 

The last equation (|65]) can be compared, after rescal- 
ing, to the dimensional equation (|8]), which we rewrite 
here for convenience after some rearrangement and sim- 
plification (e.g. 732 = 0): 



X + (2711 + 73ix^)x + {ujq + a3x'^)x 



CoVr^ 



DC 



cos ut 



2md 



(1-1) 



(66) 



The comparison of Eq. (|65|) with Eq. (|66|) results in: 



711 
731 

F = 



d^ 47VL2 

UJo_ 

_ ^711 



2md 



I ' 
fjdujQ. 



(67b) 

(67c) 
(67d) 



The last results can be used to estimate the lower bound 
of nonlinear damping due to nonlinear pretension of a 
viscoelastic string. Using Eqs. (pT|) . (|62|) . (|67|) . and 



Ah" 
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for prismatic cross-section, one has 



Pn 



2ui3 8 1 fL\ N 



3 f3 TT^ Q \h J EA Sir^Qa 



(68) 



where h denotes the dimension of the beam-string in the 
transverse direction w. 

It is possible to estimate the order of magnitude of 
Pmin in Eq. (|68|) for metals using the fact that the Young 
modulus ofbulk metals 0(10^0)^0(10^1) Pa. Also, 
the largest value of N/A that is still compatible with 
elastic behavior can be approximated by half the ulti- 
mate tensile strength, which is about 50 ^ 100 x 10^ Pa 
for most metals. For our beam-strings discussed above, 
L = 100 ^ 200 pm, h ^ 1 pm. Using these values re- 
sults in p ^ O(10~^) ^ O(10~^). For longer and wider 
beams {L = 500 pm, = 1.5 pm ) fa bricated and mea- 
sured using the same methods [90|, the lower bound 
on nonlinear damping coefficient given by Eq. (|68|) is 
Pmin ~ 0.0022 ^ 0.045, while the range of values ex- 
tracted from the experiment is 0.015 < p < 0.151 [9Q| . 
Although the elastic properties of a specific metal or alloy 
used in micro machined devices might differ significantly 
from the bulk values, they are still likely to fall inside 
the ranges defined above. Therefore, a linear viscoelastic 
process with a pure Voigt-Kelvin dissipation model can 
serve as a possible lower bound but cannot account for 
the main part of nonlinear dissipation rate found in our 
experiments. 

Unfortunately, theory describing the processes under- 
lying nonlinear damping in micromechanical beam is vir- 
tually non-existent at this moment, and no clear ten- 
dencies in the value of p were observed during the ex- 
periments. Therefore, the exact behavior of nonlinear 
damping term during beam scaling and its dependence 
on the linear Q of the structure remains elusive. Further 
experiments with wider range of micromechanical beams 
are needed to establish this behavior and to pinpoint the 
most significant mechanisms of dissipation. 
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VI. SUMMARY 

In this study, the nonhnear dynamical behavior of 
an electricahy excited micromechanical doubly clamped 
beam oscillator was investigated in vacuum. The mi- 
cromechanical beam was modeled as a Duffing-like sin- 
gle degree of freedom oscillator, nonlinearly coupled to a 
thermal bath. Using the method of multiple scales, we 
were able to construct a detailed model of slow envelope 
behavior of the system, including effective noise terms. 

It follows from the model that nonlinear damping plays 
an important role in the dynamics of the micromechan- 
ical beam oscillator. Several methods for experimental 
evaluation of the nonlinear damping contribution were 
proposed, applicable at different experimental situations. 
These methods were compared experimentally and shown 
to provide similar results. The experimental values of the 
nonlinear damping constant are non negligible for all the 
beams measured. 

Also, the slow envelope model was used to describe the 
behavior of the system close to bifurcation points in the 
presence of nonlinear damping. In the vicinity of these 
points, the dynamics of the system is significantly slowed 
down, and the phase plane motion becomes essentially 
one-dimensional. We have defined several parameters 
that govern the dynamics of the micromechanical beam 
oscillator in these conditions, and have provided simple 
approximations that can be used to estimate these pa- 
rameters from experimental data. 

The approximations developed in this study can be uti- 
lized by the experimentalist in order to estimate the in- 
accuracy of frequency response measurements of Duffing- 
like oscillators in the vicinity of bifurcation points. Ap- 
plying these results to our samples, we have found that 
thermal escape process near the bifurcation point causes 
measurement inaccuracy that is negligible. In contrast, 
the slowing down phenomenon, which is a characteristic 
of saddle-node bifurcation, can contribute a significant 



error to the measured frequency response. This error is 
non negligible even at relatively slow frequency sweeping 
rates. Similar methods can be utilized for other parame- 
ter sweeping measurements, such as excitation amplitude 
sweeping. 

As part of an effort to explain the origins of the nonlin- 
ear damping, we have formulated and analyzed a model 
of a planar, weakly nonlinear pretensioned, viscoelas- 
tic string augmented by linear Euler-Bernoulli bending, 
which incorporates a Voigt-Kelvin constitutive relation- 
ship. This model exemplifies one of the possible causes of 
non negligible nonlinear damping observed in the experi- 
ment. Based on this model, we have determined a simple 
relation connecting the maximal expected value of the 
nonlinear damping parameter, the bulk Young modulus 
of the material, and its yielding stress. However, while 
this model can serve as a lower bound, it cannot account 
for the full magnitude of the nonlinear damping mea- 
sured in the experiment. Additional experimental and 
theoretical work is required to enhance our understand- 
ing of the phenomenon of nonlinear damping in micro- 
electromechanical systems. 

In this work we have demonstrated conclusively that 
nonlinear damping in micromechanical doubly-clamped 
beam oscillator may play an important role. The meth- 
ods presented in this paper may allow a systematic study 
of nonlinear damping in micro- and nanomechanical oscil- 
lators, which may help revealing the underlying physical 
mechanisms. 
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